#include "global.h"
#include <bits/stdc++.h>
#include "solve.cpp"
//#include <eigen3/Eigen/Eigen>

int calcPropellant(void);
int calcStdHSG(double t);

int calcCmp(double T);
int findT(void);
int findTAcr(double T1,double T2);
int calcCpCv(void);
int calcGamma(void);
int findTThroatAcr(double T1,double T2);
//double calcGammaMean(double gamma);
int calcFrzCpCv(void);
int findTExitAcr(double T1,double T2);

int calcPropellant(void)
{
    cout<<"Calc Propellant."<<endl;
    compoundData tempPpl;
    //Calc atoNum
    {
        cout<<"AtoNum(1kg):"<<endl;
        //Init
        for(int j=0;j<eleNum;j++){
            tempPpl.atoNum[j]=0;
        }
        //Fuel
        for(int i=0;i<fuelNum;i++){
            for(int j=0;j<eleNum;j++){
                tempPpl.atoNum[j]+=fuel[i].atoNum[j]*fuel[i].mol;
            }
        }
        //Oxidizer
        for(int i=0;i<oxidizerNum;i++){
            for(int j=0;j<eleNum;j++){
                tempPpl.atoNum[j]+=oxidizer[i].atoNum[j]*oxidizer[i].mol;
            }
        }
        //Calc Molecular Weight
        for(int j=0;j<eleNum;j++){
            tempPpl.mass+=tempPpl.atoNum[j]*ele[j].mass;
        }
        //Scale to 1kg
        for(int j=0;j<eleNum;j++){
            propellant.atoNum[j]=tempPpl.atoNum[j]/tempPpl.mass;
            if(tempPpl.atoNum[j]!=0){
                cout<<ele[j].name<<":\t"<<propellant.atoNum[j]<<endl;
            }
        }
    }
    //Calc H
    {
        cout<<"Enthalpy=";
        //Init
        tempPpl.h=0;
        //Fuel
        for(int i=0;i<fuelNum;i++){
            tempPpl.h+=fuel[i].h*fuel[i].mol;
        }
        //Oxidizer
        for(int i=0;i<oxidizerNum;i++){
            tempPpl.h+=oxidizer[i].h*oxidizer[i].mol;
        }
        //Scale to 1kg
        propellant.h=tempPpl.h/tempPpl.mass;
        cout<<propellant.h<<"J/kg"<<endl;
    }

    cout<<"Done calc propellant."<<endl;
    return 0;
}

int calcStdHSG(double T)
{
    int j;
    for(int i=0;i<compoundNum;i++){
        for(int k=compound[i].itnNum-1;k>-1;k--){
            if(T>compound[i].tLB[k]){
                j=k;
                //cout<<"T internal="<<j<<endl;
                break;
            }
            //T internal
        }
        compound[i].stdH=R*T*(
            -pow(T,-2)/1*compound[i].a[j][0]+
            pow(T,-1)*compound[i].a[j][1]*log(T)+
            compound[i].a[j][2]+
            pow(T,1)/2*compound[i].a[j][3]+
            pow(T,2)/3*compound[i].a[j][4]+
            pow(T,3)/4*compound[i].a[j][5]+
            pow(T,4)/5*compound[i].a[j][6]+
            pow(T,-1)*compound[i].b[j][0]
            );
        compound[i].stdS=R*(
            -pow(T,-2)/2*compound[i].a[j][0]+
            -pow(T,-1)*compound[i].a[j][1]+
            log(T)*compound[i].a[j][2]+
            pow(T,1)*compound[i].a[j][3]+
            pow(T,2)/2*compound[i].a[j][4]+
            pow(T,3)/3*compound[i].a[j][5]+
            pow(T,4)/4*compound[i].a[j][6]+
            compound[i].b[j][1]
            );
        compound[i].stdMu=compound[i].stdH-compound[i].stdS*T;
        compound[i].Cp=R*(
            pow(T,-2)*compound[i].a[j][0]+
            pow(T,-1)*compound[i].a[j][1]+
            compound[i].a[j][2]+
            pow(T,1)*compound[i].a[j][3]+
            pow(T,2)*compound[i].a[j][4]+
            pow(T,3)*compound[i].a[j][5]+
            pow(T,4)*compound[i].a[j][6]
        );
    }
    return j;
}

int calcInit(void)
{
    cout<<"Calc init value."<<endl;
    for(int i=0;i<equ.N;i++){
        for(int j=0;j<eleNum;j++){
            if(compound[equ.cmpIdx[i]].atoNum[j]!=0){
                equ.c[i]=propellant.atoNum[j]/compound[equ.cmpIdx[i]].atoNum[j]/20;
                cout<<compound[equ.cmpIdx[i]].name<<":"<<equ.c[i]<<endl;
                break;
            }
        }
    }
    return 0;
}

int calcAlpha(void)
{
    //calc excess oxygen coefficient
    return 0;
}

int calcEleCmpToBeCalc(void)
{
    //Ele to be calced;
    {
        cout<<"Ele:\t";
        for(int i=0;i<eleNum;i++){
            if(propellant.atoNum[i]!=0){
                equ.ele[equ.M].name=ele[i].name;
                equ.ele[equ.M].mass=ele[i].mass;
                equ.eleIdx[equ.M]=i;
                cout<<equ.ele[equ.M].name<<" ";
                equ.M+=1;
            }
        }
        cout<<endl;
    }
    //Compound to be cacled.
    {
        cout<<"Compound:\t";
        for(int i=0;i<compoundNum;i++){
            compound[i].toBeCalc=true;
            for(int j=0;j<eleNum;j++){
                if(compound[i].atoNum[j]!=0){
                    if(propellant.atoNum[j]==0){
                        compound[i].toBeCalc=false;
                        break;
                    }
                }
            }
            if(compound[i].toBeCalc)
                cout<<compound[i].name<<" ";
        }
        cout<<endl;
        equ.L=equ.N=0;
        for(int i=0;i<compoundNum;i++){
            if(compound[i].toBeCalc==true){
                if(compound[i].phase==0){
                    equ.gasIdx[equ.N-equ.L]=i;
                }
                else{
                    equ.cdsIdx[equ.L]=i;
                    equ.L+=1;
                }
                equ.N+=1;
            }
        }
        for(int i=0;i<equ.L;i++){
            equ.cmpIdx[i]=equ.cdsIdx[i];
            //cout<<equ.cmpIdx[i]<<" ";
        }
        for(int i=equ.L;i<equ.N;i++){
            equ.cmpIdx[i]=equ.gasIdx[i-equ.L];
            //cout<<equ.cmpIdx[i]<<" ";
        }
        cout<<endl;
    }
    return 0;
}

int calcMat(void)
{
    cout<<"L="<<equ.L<<"\tN="<<equ.N<<"\tM="<<equ.M<<endl;
    //calc A[N][L]
    {
        cout<<"A="<<endl;
        for(int i=0;i<equ.M;i++){
            cout<<"\t";
            for(int j=0;j<equ.N;j++){
                equ.A[i][j]=compound[equ.cmpIdx[j]].atoNum[equ.eleIdx[i]];
                cout<<equ.A[i][j]<<" ";
            }
            cout<<endl;
        }
    }
    //calc R[M][M]
    {
        cout<<"R="<<endl;
        for(int k=0;k<equ.M;k++){
            cout<<"\t";
            for(int i=0;i<equ.M;i++){
                equ.R[k][i]=0;
                for(int j=equ.L;j<equ.N;j++){
                    equ.R[k][i]+=equ.A[k][j]*equ.A[i][j]*equ.c[j];
                }
                cout<<equ.R[k][i]<<"\t";
            }
            cout<<endl;
        }
    }
    //calc a[M]
    {
        cout<<"a="<<endl<<"\t";
        for(int k=0;k<equ.M;k++){
            equ.a[k]=0;
            for(int j=equ.L;j<equ.N;j++){
                equ.a[k]+=equ.A[k][j]*equ.c[j];
            }
            cout<<equ.a[k]<<"\t";
        }
        cout<<endl;
    }
    //calc mu[N}]
    {
        cout<<"mu="<<endl;
        calcStdHSG(equ.T);
        for(int j=0;j<equ.N;j++){
            equ.mu[j]=compound[equ.cmpIdx[j]].stdMu;
            cout<<"\t"<<equ.mu[j];
        }
        cout<<endl;
    }
    //calc y[N]
    {
        cout<<"y="<<endl;
        for(int j=0;j<equ.N;j++){
            equ.y[j]=-equ.mu[j]/(R*equ.T);
            cout<<"\t"<<equ.y[j];
        }
        cout<<endl;
    }
    //calc cG
    {
        cout<<"c_g="<<endl;
        equ.cG=0;
        for(int j=equ.L;j<equ.N;j++){
            equ.cG+=equ.c[j];
        }
        cout<<"\t"<<equ.cG<<endl;
    }
    //calc Phi
    {
        cout<<"Phi="<<endl;
        for(int j=0;j<equ.L;j++){
            equ.Phi[j]=-equ.y[j]*equ.c[j];
            cout<<"\t"<<equ.Phi[j];
        }
        for(int j=equ.L;j<equ.N;j++){
            equ.Phi[j]=(-equ.y[j]+log(equ.c[j])+log(equ.p/pAtm)-log(equ.cG))*equ.c[j];
            cout<<"\t"<<j<<":"<<equ.Phi[j];
        }
        cout<<endl;
    }
    //calc vecB
    {
        for(int k=0;k<equ.M;k++){
            equ.vecB[k]=propellant.atoNum[equ.eleIdx[k]];
            //cout<<equ.vecB[k]<<"\t";
            //equ.vecB[k]=0;
            //for(int j=0;j<equ.N;j++){
            //    equ.vecB[k]+=equ.A[k][j]*equ.c[j];
            //}
            for(int j=equ.L;j<equ.N;j++){
                equ.vecB[k]+=equ.A[k][j]*equ.Phi[j];
            }
            //cout<<endl;
        }
        equ.vecB[equ.M]=0;
        for(int j=equ.L;j<equ.N;j++){
            equ.vecB[equ.M]+=equ.Phi[j];
        }
        for(int k=equ.M+1;k<equ.M+equ.L+1;k++){
            equ.vecB[k]=-equ.y[k-equ.M-1];
        }
        //cout
        cout<<"vecB="<<endl;
        for(int k=0;k<equ.M+equ.L+1;k++){
            cout<<"\t"<<equ.vecB[k];
        }
        cout<<endl;
    }
    //calc matA
    {
        for(int i=0;i<equ.M;i++){
            for(int j=0;j<equ.M;j++){
                equ.matA[i][j]=equ.R[i][j];
            }
        }
        for(int i=0;i<equ.M;i++){
            equ.matA[i][equ.M]=equ.matA[equ.M][i]=equ.a[i];
        }
        for(int i=0;i<equ.M;i++){
            for(int j=equ.M+1;j<equ.M+equ.L+1;j++){
                equ.matA[i][j]=equ.matA[j][i]=equ.A[i][j-equ.M-1];
            }
        }
        //cout
        cout<<"matA="<<endl;
        for(int i=0;i<equ.M+equ.L+1;i++){
            for(int j=0;j<equ.M+equ.L+1;j++){
                cout<<"\t"<<equ.matA[i][j];
            }
            cout<<endl;
        }
    }
    return 0;
}

int calcHF(void){
    equ.hF=0;
    double temp=0;
    for(int j=0;j<equ.N;j++){
        equ.hF+=compound[equ.cmpIdx[j]].stdH*equ.X[j];
    }
    for(int j=equ.L;j<equ.N;j++){
        temp+=equ.X[j];
    }
    //cout<<"h_f="<<endl<<"\t"<<equ.hF<<endl;
    equ.molWeight=1/temp;
    //cout<<"Molecular weight="<<"\t"<<equ.molWeight<<"kg/mol"<<endl;
    return 0;
}

int calcCmp(double T){
    equ.T=T;
    calcMat();
    solveEqu();
    for(int i=0;i<15;i++){
        updateInit();calcMat();solveEqu();
    }
    calcHF();
    return 0;
}

int findT(void)
{
    double hF1,hF2;
    double T1,T2;
    T1=2000;
    calcCmp(T1);
    hF1=equ.hF-propellant.h;
    for(double T=2100;T<=4500;T+=100){
        calcCmp(T);T2=T;
        hF2=equ.hF-propellant.h;
        if((hF1>0&&hF2<=0)||(hF1<=0&&hF2>0)){
            //TF=T1+(-hF1)/(hF2-hF1)*(T2-T1);
            //cout<<"T_f="<<endl<<"\t"<<TF<<endl;
            findTAcr(T1,T2);
            return 0;
        }
        hF1=hF2;
        T1=T2;
    }
    return 0;
}

int findTThroat(void)
{
    double hF1,hF2;
    double T1,T2;
    T1=2000;
    calcCmp(T1);
    //calc properties
    calcCpCv();
    calcGamma();
    hF1=equ.hF-propellant.h+equ.Gamma*R/equ.molWeight*T1/2;
    for(double T=2100;T<=4500;T+=100){
        calcCmp(T);T2=T;
        hF2=equ.hF-propellant.h+equ.Gamma*R/equ.molWeight*T/2;
        if((hF1>0&&hF2<=0)||(hF1<=0&&hF2>0)){
            //TF=T1+(-hF1)/(hF2-hF1)*(T2-T1);
            //cout<<"T_f="<<endl<<"\t"<<TF<<endl;
            findTThroatAcr(T1,T2);
            return 0;
        }
        hF1=hF2;
        T1=T2;
    }
    return 0;
}

int findTExit(void)
{
    double hF1,hF2;
    double T1,T2;
    double k=nozzleThroat.Gamma;
    T1=1500;
    calcCmp(T1);
    double piRatio=nozzleExit.p/nozzleThroat.p;
    nozzleExit.lambda=1-pow(piRatio,(k-1)/k)*(1-(k-1)/(k+1));
    nozzleExit.lambda/=(k-1)/(k+1);
    nozzleExit.lambda=sqrt(nozzleExit.lambda);
    //calc properties
    calcCpCv();
    calcGamma();
    hF1=equ.hF-propellant.h+equ.Gamma*R/equ.molWeight*T1/2*pow(nozzleExit.lambda,2)*2/(k+1-(k-1)*pow(nozzleExit.lambda,2));
    for(double T=2100;T<=4500;T+=100){
        calcCmp(T);T2=T;
        double piRatio=nozzleExit.p/nozzleThroat.p;
        nozzleExit.lambda=1-pow(piRatio,(k-1)/k)*(1-(k-1)/(k+1));
        nozzleExit.lambda/=(k-1)/(k+1);
        nozzleExit.lambda=sqrt(nozzleExit.lambda);
        hF2=equ.hF-propellant.h+equ.Gamma*R/equ.molWeight*T/2*pow(nozzleExit.lambda,2)*2/(k+1-(k-1)*pow(nozzleExit.lambda,2));
        if((hF1>0&&hF2<=0)||(hF1<=0&&hF2>0)){
            //TF=T1+(-hF1)/(hF2-hF1)*(T2-T1);
            //cout<<"T_f="<<endl<<"\t"<<TF<<endl;
            findTExitAcr(T1,T2);
            return 0;
        }
        hF1=hF2;
        T1=T2;
    }
    return 0;
}

int findTExitAcr(double T1,double T2)
{
    if(T2-T1<eps2){
        TF=(T2+T1)/2;
        cout<<"T_f="<<endl<<"\t"<<TF<<endl;
        return 0;
    }
    double hF1,hF2;
    calcCmp(T1);
    double k=nozzleThroat.Gamma;
    double piRatio=nozzleExit.p/nozzleThroat.p;
    nozzleExit.lambda=1-pow(piRatio,(k-1)/k)*(1-(k-1)/(k+1));
    nozzleExit.lambda/=(k-1)/(k+1);
    nozzleExit.lambda=sqrt(nozzleExit.lambda);
    hF1=equ.hF-propellant.h+equ.Gamma*R/equ.molWeight*T1/2*pow(nozzleExit.lambda,2)*2/(k+1-(k-1)*pow(nozzleExit.lambda,2));
    calcCmp(T2);
    nozzleExit.lambda=1-pow(piRatio,(k-1)/k)*(1-(k-1)/(k+1));
    nozzleExit.lambda/=(k-1)/(k+1);
    nozzleExit.lambda=sqrt(nozzleExit.lambda);
    hF2=equ.hF-propellant.h+equ.Gamma*R/equ.molWeight*T2/2*pow(nozzleExit.lambda,2)*2/(k+1-(k-1)*pow(nozzleExit.lambda,2));
    TF=T1+(-hF1)/(hF2-hF1)*(T2-T1);
    if(abs(TF-T1)>abs(TF-T2)){
        findTExitAcr((T1+TF)/2,T2);
    }
    else{
        findTExitAcr(T1,(T2+TF)/2);
    }
    return 0;
}

int findTThroatAcr(double T1,double T2)
{
    if(T2-T1<eps2){
        TF=(T2+T1)/2;
        cout<<"T_f="<<endl<<"\t"<<TF<<endl;
        return 0;
    }
    double hF1,hF2;
    calcCmp(T1);
    hF1=equ.hF-propellant.h+equ.Gamma*R/equ.molWeight*T1/2;
    calcCmp(T2);
    hF2=equ.hF-propellant.h+equ.Gamma*R/equ.molWeight*T2/2;
    TF=T1+(-hF1)/(hF2-hF1)*(T2-T1);
    if(abs(TF-T1)>abs(TF-T2)){
        findTThroatAcr((T1+TF)/2,T2);
    }
    else{
        findTThroatAcr(T1,(T2+TF)/2);
    }
    return 0;
}

int findTAcr(double T1,double T2){
    if(T2-T1<eps){
        TF=(T2+T1)/2;
        cout<<"T_f="<<endl<<"\t"<<TF<<endl;
        return 0;
    }
    double hF1,hF2;
    calcCmp(T1);
    hF1=equ.hF-propellant.h;
    calcCmp(T2);
    hF2=equ.hF-propellant.h;
    TF=T1+(-hF1)/(hF2-hF1)*(T2-T1);
    if(abs(TF-T1)>abs(TF-T2)){
        findTAcr((T1+TF)/2,T2);
    }
    else{
        findTAcr(T1,(T2+TF)/2);
    }
    return 0;
}

int calcCpCv(){
    double T0,T1,T2,p1,p2,p0;
    T0=TF;
    T1=TF+1;
    T2=TF-1;
    p1=equ.p+100;
    p2=equ.p-100;
    p0=equ.p;
    //calc Cp
    double h1,h2,M1,M2,M3,M4;
    equ.p=p0;
    calcCmp(T1);
    h1=equ.hF;
    M1=equ.molWeight;
    calcCmp(T2);
    h2=equ.hF;
    M2=equ.molWeight;
    equ.cP=(h2-h1)/(T2-T1);
    //calc Cv
    equ.p=p1;
    calcCmp(T0);
    M3=equ.molWeight;
    equ.p=p2;
    calcCmp(T0);
    M4=equ.molWeight;
    double A;
    A=1-(log(M1)-log(M2))/(log(T1)-log(T2));
    A*=A;
    A/=(1+(log(M3)-log(M4))/(log(p1)-log(p2)));
    equ.cV=equ.cP-A*R/equ.molWeight;
    return 0;
}

int calcGamma(void)
{
    equ.Gamma=equ.cP/equ.cV;
    return 0;
}

int calcChamber(void){
    //redirect
    ofstream fout("calc.out");
    streambuf *oldcout;
    oldcout = cout.rdbuf(fout.rdbuf());
    //calc
    equ.p=chamber.p;
    calcInit();
    findT();
    //calc properties
    calcCpCv();
    calcGamma();
    //calc chamber fractions
    calcCmp(TF);
    solveEqu();
    //velocity
    equ.velocity=0;
    //redirect
    cout.rdbuf(oldcout);
    //trans
    chamber.T=TF;
    chamber.molWeight=equ.molWeight;
    chamber.Gamma=equ.Gamma;
    return 0;
}

int calcNozzleThroat(void)
{
    //redirect
    ofstream fout("calc.out");
    streambuf *oldcout;
    oldcout = cout.rdbuf(fout.rdbuf());
    //calc
    double k=equ.Gamma;
    equ.p=chamber.p*pow(1+(k-1)/2,-k/(k-1));
    TF*=1-(k-1)/(k+1);
    //find T
    findTThroat();
    //calc properties
    calcCpCv();
    calcGamma();
    //calc nozzle throat fractions
    calcCmp(TF);
    solveEqu();
    //velocity
    equ.velocity=sqrt(equ.Gamma*R/equ.molWeight*TF);
    //redirect
    cout.rdbuf(oldcout);
    //trans
    nozzleThroat.p=equ.p;
    nozzleThroat.T=TF;
    nozzleThroat.molWeight=equ.molWeight;
    nozzleThroat.Gamma=equ.Gamma;
    nozzleThroat.h=equ.hF;
    nozzleThroat.velocity=equ.velocity;
    return 0;
}

int calcNozzleExitEqm(void)
{
    //redirect
    ofstream fout("calc.out");
    streambuf *oldcout;
    oldcout = cout.rdbuf(fout.rdbuf());
    //calc
    double k=equ.Gamma;
    equ.p=nozzleExit.p;//*pow(1+(k-1)/2,-k/(k-1));
    double lambda;
    double piRatio=nozzleExit.p/nozzleThroat.p;
    lambda=1-pow(piRatio,(k-1)/k)*(1-(k-1)/(k+1));
    lambda/=(k-1)/(k+1);
    lambda=sqrt(lambda);
    TF/=1-(k-1)/(k+1);
    TF*=1-(k-1)/(k+1)*pow(lambda,2);
    //find T
    findTExit();
    //calc properties
    calcCpCv();
    calcGamma();
    //calc nozzle throat fractions
    calcCmp(TF);
    solveEqu();
    calcCpCv();
    calcGamma();
    piRatio=nozzleExit.p/nozzleThroat.p;
    lambda=1-pow(piRatio,(k-1)/k)*(1-(k-1)/(k+1));
    lambda/=(k-1)/(k+1);
    lambda=sqrt(lambda);
    //velocity
    //equ.velocity=sqrt(equ.Gamma*R/equ.molWeight*TF);
    equ.velocity=sqrt(2*(nozzleThroat.h-equ.hF)+pow(nozzleThroat.velocity,2));
    //redirect
    cout.rdbuf(oldcout);
    //trans
    nozzleExit.T=TF;
    nozzleExit.molWeight=equ.molWeight;
    nozzleExit.Gamma=equ.Gamma;
    nozzleExit.h=equ.hF;
    nozzleExit.velocity=equ.velocity;
    return 0;
}

int calcNozzleExit(void)
{
    //calc
    double k=nozzleThroat.Gamma;
    TF=nozzleThroat.T*pow((nozzleThroat.p/nozzleExit.p),(1-k)/k);
    calcStdHSG(TF);
    calcFrzCpCv();
    calcGamma();
    //update k
    nozzleExit.T=TF;
    k=(k+equ.Gamma)/2;
    TF=nozzleThroat.T*pow((nozzleThroat.p/nozzleExit.p),(1-k)/k);
    double meanGamma;
    for(int i=0;i<1;i++){
        meanGamma=log(nozzleExit.p/chamber.p)/log(nozzleExit.p/chamber.p*nozzleThroat.molWeight/chamber.molWeight*chamber.T/TF);
        TF=chamber.T*pow(nozzleExit.p/chamber.p,(meanGamma-1)/meanGamma);
    }
    calcStdHSG(TF);
    /*
    calcStdHSG(TF);
    calcFrzCpCv();
    calcGamma();
    */
    calcHF();
    equ.velocity=sqrt(2*(nozzleThroat.h-equ.hF)+pow(nozzleThroat.velocity,2));
    //equ.velocity=
    /*
    calcStdHSG(TF);
    calcFrzCpCv();
    calcGamma();
    */
    equ.p=nozzleExit.p;
    nozzleExit.velocity=equ.velocity;
    nozzleExit.molWeight=nozzleThroat.molWeight;
    nozzleExit.velocity=equ.velocity;
    nozzleExit.T=TF;
    return 0;
}

int calcFrzCpCv(void)
{
    equ.cP=0;
    double temp=0;
    for(int j=0;j<equ.N;j++){
        equ.cP+=compound[equ.cmpIdx[j]].Cp*equ.X[j];
    }
    equ.cV=equ.cP-R/equ.molWeight;
    //cout<<"Cp="<<endl<<"\t"<<equ.cP<<endl;
    return 0;
}

/*
double calcGammaMean(double gamma)
{
    cout<<gamma<<endl;
    double pe=nozzleExit.p;
    double pt=nozzleThroat.p;
    double Tt=nozzleThroat.T;
    double Te=Tt*pow((pe/pt),(gamma-1)/gamma);
    double newgamma=1;
    newgamma/=1-log(Tt/Te)/log(pt/pe);
    if(abs(gamma-newgamma)<eps){
        cout<<newgamma<<endl;
        return newgamma;
    }
    else return calcGammaMean(newgamma);
}*/

int calcPerf(void)
{
    vcmImpulse=nozzleExit.velocity+R/nozzleExit.molWeight*nozzleExit.T/nozzleExit.velocity;
    double GGamma,k;
    k=chamber.Gamma;
    GGamma=sqrt(k)*pow(2/(k+1),(k+1)/(2*(k-1)));
    charVel=sqrt(R/nozzleExit.molWeight*chamber.T)/GGamma;
    return 0;
}
